
#include<iostream>
using namespace std;

#include<cmath>
#include<ctime>

#include"library.h"

int main(){

	srand(time(NULL));
	int i,j,k;
	int n=3;
	double y;

	double * a = new double [n*n];
	double *a_diag = new double [n*n];

	for(i=0;i<n;i++){
		for(j=i;j<n;j++){
			a[ i*n + j ] = (i+2)*(j+2)+3;//(double)rand()/RAND_MAX;
			a[ j*n + i ] = a[ i*n + j];
			a_diag[ i*n + j] = a[ i*n + j];// a_diag = a
			a_diag[ j*n + i] = a[ j*n + i];
		}
	}
	cout<<"a:"<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			a_diag[ i*n+j ] = a[ i*n+j ];
			cout<<a_diag[ i*n+j ]<<",";
		}
		cout<<endl;
	}

	double *z = new double [n];

	jacobi(n, a_diag, z);
	cout<<"jacobi eigenvalues: "; for(int i=0;i<n;i++)cout<<z[i]<<","; cout<<endl;

	cout<<"A v: "<<endl;
	for(i=0;i<n;i++){// ith column of ( A (a_diag)^top )
		cout<<" A v"<<i<<" = ";
		for(j=0;j<n;j++){
			y=0;
			for(k=0;k<n;k++) // A \vec{v} 
				y += a[j*n+k] * a_diag[i*n+k];
			cout<<y<<",";
		}
		cout<<endl;
	}
	cout<<" lambda * v : "<<endl;
	for(i=0;i<n;i++){
		cout<<" lambda"<<i<<" * v"<<i<<" = ";
		for(j=0;j<n;j++)cout<< z[i]* a_diag[ i*n + j]<<",";
		cout<<endl;
	}

	delete [] z,a,a_diag;
	return 0;
}
